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Abstract 

We consider a stochastic model of the two-dimensional chemostat as a 
diffusion process for the concentration of substrate and the concentration of 
biomass. The model allows for the washout phenomenon: the disappear- 
ance of the biomass inside the chemostat. We establish the Fokker-Planck 
associated with this diffusion process, in particular we describe the boundary 
conditions that modelize the washout. We propose an adapted finite difference 
scheme for the approximation of the solution of the Fokker-Planck equation. 
Keywordsichemostat, stochastic differential equation, Fokker-Planck equa- 
tion, finite difference scheme 

1 Introduction 

Many biotechnological processes are modelized with the help of ordinary differential 
equations (ODE). For example, the dynamic for a single species/single substrate 
chemostat is classically modelized as [12]: 

s{t) = -k yu(s(t)) b{t) + D - s{t)) , (la) 
b{t) = {^^{s{t))-D}b{t) (lb) 

where b{t) and s{t) are the concentrations of biomass and substrate at time t inside 
the chemostat. The parameters are the dilution rate D, the input substrate concen- 
tration Sin, and the stoichiometric coefficient k. The specific growth function fi{s) 
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could be of the Monod (non-inhibitory) type: 

where /i^ax is the maximum growth rate and ks is the half-saturation; it could also 
be of the Haldane (inhibitory) type: 

K^) = 1 o ; • (3) 

^ ' k, + s + s^/a ^ ' 

As pointed out in [2], the system (1) is simple and applicable to many situations, 
it can be seen as a limit model of a stochastic birth and death process in high 
population size asymptotic. Hence (1) can give account for the mean behavior of 
the underlying stochastic process but it cannot give account for its the variance. 
Moreover (1) fails to propose a realistic representation of the chemostat in small 
population scenario, that is in cases close to the washout (corresponding to the 
disappearance of the biomass, i.e. hif) = 0). 

We present the stochastic model in Section 2 and derive the associated Fokker- 
Planck equation in Section 3. A finite difference scheme approximation is detailled 
in Section 4 and some numerical tests are presented in Section 5. 



2 The stochastic chemostat model 

Consider the stochastic process Xt = {X^,Xl) = {St,Bt) solution of: 

d5i = { - kfi{St) Bt + D (s„, - St)}dt + ci ^/s'tdWl , (4a) 
dBt = MSt) -D}Btdt + C2^t dW^ , (4b) 

where Bt and St are the concentrations of biomass and substrate at time t; and 
Wl are independent scalar standard Brownian motions; ci > and C2 > are the 
noise intensities; Wl and W"^ are independent scalar standard Wiener processes. We 
suppose that 5*0 > and > so that S'j > and -Bt > for alH > 0. 

The precise analysis of the behavior of the solution of (4) will be addressed in 
a forthcoming work [3]. Still we can describe it simply with some highlights about 
the classic Cox-lngersoU-Ross model. Consider the one-dimensional SDE: 

dit = {a + hit)dt + a^tdWt, = > . (5) 

with a > 0, 6 G M, cr > 0. According to [10, Prop. 6.2.4], for all xq > 0, is a 
continuous process taking values in M"*", and let r = mi{t > 0, = 0}, then: 

(z) If a > cr^/2, then r = oo Pa;-a.s.; 
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(a) if < a < and 6 < then r < oo F^-a.s.; 
{m) if < a < aV2 and 6 > then P^(r < oo) G (0, 1). 

In the first case, S,t never reaches 0. In the second case S,t a-s. reaches the state 0, in 
the third case it may reach 0. If a = then the state is absorbing. 

In case of the System (4), it is clear that i? = is an absorbing state for (4b), and 
when 5 = 0, (4a) reduces to the substrate dynamics conditionally of the washout, 
namely: 

dSl = D (Si„ - S;) dt + ci dWl (6) 
hence the solution of this SDE will stay on the half-line [0, oo) and: 
(z) if D Sin > ^ then St never reaches 0; 

(a) a D Sin < ^ then St reaches in finite time and is reflected. 

Note that, as Ci is "small", condition (z) is more realistic than condition (ii): indeed, 
with a continuous input Si„ , there is no reason for the substrate concentration in 
the chemostat to vanish. 

Simulation schemes for (4) should respect the previous properties, an adequate 
choice is: 

St+5 =[St + {-k fi{St) Bt + D (sin -St)]5 + ci ^t^5w]] _^ , (7a) 
Bt+5 = [Bt + {^^{St) -D}Bt5 + C2^t^5w^]^, (7b) 

where {Wj^^jigN and {w1^}ii^n are i.i.d. A^(0, 1) random variables, also independent 
from Xq. Note that i?f = is absorbing for (7b). 

Notations 2.1 Let x = (xi,X2) = (s, 6) G = [0,oo)^ and 

h{x) = fi{s,b) = -kfi{s)b + D (s,„ - s) , ai{x) = (Ti(s,6) = ai{s) = ci^/s, 
f2{x) = f2{s,b) = [/i(s) - D]b, a2{x) = 0-2(5,6) = 0-2(6) = 02^6, 

so that (4) reads: 

dXt = f{Xt)dt + a{Xt)dWt 

^^th fix) ={iU)' -(^) = ( ^l) )andWt={^^). 

Let dRl = Ti U 12 with Ti = {{s,b) G [0,oo)^ 6 = 0} and 12 = {(s,6) G 
[0,oo)2; s = 0}. 
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3 The Fokker-Planck equation 

Let 7!'t{dx) = 7!'t{ds, db) be the distribution law of of Xt = {St, Bt): 

iVt{A,B)=F{SteA, BteB) 
for all Borel sets A, B of [0, oo). According to [11], iitidx) of Xt can be decomposed 



as: 



Tit^dx) = nt{ds X db) = So{db) qt{s) ds + pt{s, b) ds db (8) 

indeed the diffusion process "lives" in but never reaches T2 so the distribution law 
features only a "regular" component b) that only charges and a "degenerate" 
component qt{s) that only charges Fi. 

As Tit is a probability distribution we get the normalization property: 



00 fOO POO 



qt{s)ds+ / / pt{s,b)dsdb = 1 
Jo Jo 



and the washout probability at time t is: 



PCO POO POO 

F{Bt = 0)= qt{s)ds = l- / pt{s,b)dsdb. 
Jo Jo Jo 



The Fokker-Planck equation in a weak form is: 



d_ 
dt 



7it{ds,db) (f){s,b) 



[[ ntids,db)C<P{s,b) (9) 
J JrI 



for all test functions 0, where C is the infinitesimal generator defined by: 

L(j){x) = C(f){s,b) 

df ^ 1 2 

= Ms, b) b) + Us, b) ct>',{s, b) + §s b) + f b) . (lo) 

Using the decomposition (8), the Fokker-Planck equation (9) reads: 
^1 r qt{s)4>{s,0)ds+ l!vt{s,b)<P{s,b)dsdb] = 

qt{s)C(t){s,Q)ds+ [[ ptis,b)C(l)is,b)dsdb (11) 
J Jm.1 
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Lemma 3.1 For all functions (f) G H^^{M.1_) (i.e. (f) e H^(U.1.) and 0|r2 = 0^ and 
t > 



„2 poo 
^2 



(pt,£0)=/ £>t(x)0(a;)dx+^ / 0) 0) ds 







•where C* is the adjoint operator: 

-[^(x) h{x)l - [^(x) + f [^(x) + f [^(x) 6];'. . 

Proof By definition of £: 

(pt, £0) = 4^ pt(x) £0(x) dx = 4^ pt(x) /i(x) 0;(x) dx + ^2 pt{x) /2(x) (p[{x) dx 

+ f /r2 ^ (3^) da; + f /ig2^ (x) 6 0'f;2 (a;) dx 

we consider separately tliese four last terms. 

From Green's formula [1]: u'^^vdx = — J^2 uv'^_dx + Jq^2 uvUidS^ where 

Ui is the ith component of the outward unit normal n, i.e. ni{x) = on Fi and —1 
on F2 and n2(x) = —1 on Fi and on F2. So we get: 

/r^ Pt{x) /i(x) 0;(x) dx = - /^^ [pt{x) /i(x)]', 0(x) dx + /gjg^ pt{x) /i(x) 0(x) ni(x) d5^ 

= - /k2 </>(a;) dx - pt(x) /i(x) 0(x) d5^. 

= - !w^_^[Pt{x) /i(x)]; </)(x) dx . (as = on F2) 

For the second term: 

/ir2 Pt{x) /2(x) (j)[{x) dx = - 4^ [pt(x) /2(a;)]'f, 0(x) dx + /^^a /2(x) 0(x) n2(x) d>S^ 

= - 4^ [Pt(a;) /2(a;)]'f, 0(x) dx - /^^ pt(x) /2(x) 0(x) d>S^ 
= - /jj2^ [Pt{.x) /2(a;)]'f, 0(x) dx . (as /2 = on Fi) 

For the third term: 

/jg2 Pt {x) s (j)'^2 (x) dx = - J^2^ [pt (x) s]; 0; (x) dx + Jg^2^ pt (x) s 0', (x) ni (x) dS^ 
= - bt(a;) s]; 0',(x) dx - /^^ pt{x) s 0',(x) d5^ 
= - J^2^ [Ptix) s]'g 0s(x) dx (as s = on F2) 

= f^2^ [pt{x) s]'^2 (p{x) dx - /gjj2 [pt{x) s]i 0(x) ni(x) d5^ 
= /m:^ bt(a;) 0(a^) da; + [pt{x) s]'^ 0(x) d5^ 
= /m2 bt(a;) s]"2 0(x) dx . (as = on F2) 
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For the fourth term: 

Jul Pti^) ^ ^2 (^) dx = - /j^2 [pt{x) 6]; dx + /g^2 6 ^'^(x) n2{x) dS^ 
= - /r2 lPt(^) (P'bi^) dx - /p^ (x) 6 0;(x) 
= - J^2_ [Pt{x) b]'^ (f)'f,{x) dx (as 6 = on Ti) 

= /k2 [Pt(x) 0(x) dx - /g^2 [pt(x) 6]; 0(x) (x) d5^ 
= /k2 [Pt{x) b]'l2 0(x) dx + [pt(x) 6]';, 0(x) d5^ . 

Summing up these identities leads to: 
finally 

/pjpt(x)6]';,0(x) dS^ = /p^ {[Pt(x)]b& + Pt(x)} 0(x) d5^ 

= /r^Pt(x) 0(x) d5^ 
= /~Pt(s,O)0(s,O) ds 

proves the lemma. □ 

According to Lemma 3.1, (11) becomes: 

^1 r qt{s)^{s,0)ds+ f[pt{s,b)(P{s,b)dsdb] = 

/•CxD /• /» 

= / gt(s)£0(s,O)ds+ // Cpt{s,b)(l){s,b)dsdb 

Jo J JR\ 

2 rco 

Let (p{s,b) = ip{s)^{b) with ^(0) = 1, ^/'(fe) = for 6 > £ and ^'(0) = ip"{0) = 0, 
after letting e — )■ 0, the previous equation leads to: 

^ /"oo poo ^2 poo 

dt 

where 



/ qt{s) (f{s) ds = qt{s)gif{s)ds + ^ / pt(s, 0) v?(s) ds (13) 
Jo Jo ^ Jo 



Q^{s) = D (s,„ - s) ^'{s) + f s <f"{s) (14) 

is the infinitesimal generator of the diffusion St in washout mode, i.e. of the SDE 
(6). As (13) is valid for all test functions if, we get the following equation for qt{s): 

g-^qtis) = g*qt{s) + ^Pt{s,0), Vt >0, [0,oo) (15a) 
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the equation for pt{s, v) is 

^^pt{s, v) = C*pt{s, v), Vt > , (s, v) G [0, oo)2 (15b) 

The initial condition for (15a) and (15b) are: 

Qti-s) = p,{s) , pt{s,v) = p{s,b) . (15c) 

where p^{s) ds So{db) + p{s, h) ds dh is the distribution law of Xq = {Sq, Bq). 
The operators are: 

GMs) = -D [(s,„ - .) ^(s)] ' + f cfis)] " , (16) 
£>(., v) = - [h{s, h) 0(5, h)] [ - [Ms, b) <P{s, b)] [ 

+ §[s<P{s,b)]':, + §[b<p{s,b)]'l, (17) 

Finally, the Fokker-Planck equation is a system of PDE's: (15b) for pt{s,v) and 
(15a) for qt{s), the first one is autonomous, and its solution appears as an input for 
the second PDE. 



4 Approximation 

Many finite difference schemes and finite element schemes are adapted to space 
discretization of the system (15). Here we use the specific finite difference scheme 
proposed in [8]. This classical scheme presents nice numerical properties and it also 
can be interpreted as an approximation of the solution of (4) by a pure jump Markov 
process on a finite discretization grid, the resulting system in discrete-space and 
continuous-time is the exact Fokker-Planck equation (forward Kolmogorov equation) 
associated with this pure jump process. The infinitesimal generator C of the SDK 
(4) is given by (10), this operator fully characterizes the distribution law of the 
process Xt = {St,Bt), indeed the set of equations (15) is totally determined by the 
operator £ as ^ is only the restriction of £ to 

The finite difference scheme is detailed in A, it leads to the following approxi- 
mation of the infinitesimal generator: 

C(f){x) ~ Ch(t)ix) = ^ Chix,y) (j){y) 
y&Gh 

for X G Gh where: 

{ki hi, A;2 /i2) ; ki = 0,..., Ni, i = 1,2} , 
{ki hi,k2h2); h = 1, Ni - 1, i = 1,2} , 
{kihi,0); ki = 0,...,Ni}, 



Gh = {x = 
Gh = {x = 
Gl {x = 
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^inax ' 



O 



subsrtate concentration 




Figure 1: Discretized domain Gh- 

are the grid version of M^, and Fi respectively, see Figure 1. 
For the interior points x G Gh the finite difference scheme is: 

^ l/lWI _ l/2(x)| _ a'iix) _ alix) 

hi h2 h'j h'j ' 

_ /±(x) aUx) _ 1 9 

~ /li ' ' — -L'^ ' 

= otherwise. 

For the boundary points x E Gh \ Gh the finite difference schemes are detailed in 
B. They correspond to the Figure 1: for s = s^^^ oi b = b^^^, we must impose 
reflecting conditions, for for s = or 6 = 0, the boundary conditions are natural, 
they derive from the value of the coefficients. Indeed, when 6 = 0, then /2 = = 
and the jump process stays on the boundary "6 = 0" (it cannot jump to 6 = /i2 or 
to 6 = —h2). When s = 0, then fi = D S;,-, and (Ti = 0, so the jump process can only 
jump to s = hi. 

We obtain a matrix Ch = ['C.h{x,y)]x,y£Gh which is the infinitesimal generator of 
a pure jump Markov process {X^)t>o in continuous time and discrete state space Gh- 
Starting from a point x of the grid, the process Xj^ stays there during a time exponen- 
tially distributed with parameter —Ch{x, x) then it jumps to a point y with probabil- 
ity Ch{x,y)/{-Ch{x,x)) for all y G Qh, and Ch is a Q-matrix as X^jyee^ '^h{x,y) = 0. 
Then the following Kolmogorov forward equation: 

^^p1{x) = CIpHx) (18) 
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gives the evolution of the distribution law of X/*, p^{x) = F{X^ = x), x G G^. 

It is important to note that this approach gives an approximation of the coupled 
system of PDEs (15): {Pt{x))x&Gh\Gl is an approximation of 6))(s,b)g(o,oo)2 and 
{Pt{^))x€Gl is an approximation of (gt(s))^6(o,oo)- 

For the time-discretization we use the implicit Euler scheme approximation: 

^ = Cf^Pt+six) 

that is: 

{I - S ^h) Pt+six) =Pt{x). 

5 Numerical results 

5.1 Comparison 

Many works [7] propose the following structure for the diffusion coefficients: 

dSt = {-k fi{St) Bt + D - St) } dt + ci St dWl , (19a) 
dBt = IfiiSt) Bt-DBt}dt + C2 Bt dW^ . (19b) 

It is slightly different from (4). In large population size, these two models are rather 
equivalent; they differ drastically in the washout regime. 

In this test we use the Monod growth rate function (2) and the parameters: 
k = 10, Si„ = 1.3 (mg/1), D = 0.4 (1/h), fi^,^ = 3 (1/h) , A;, = 6 (mg/1). The initial 
law is {So,Bo) ~ A/'(0.45, 10"^) (g) A^(0.01, 10"^). The discretization parameters are 
Smax = 2, 6,„ax = 0.06, 5 = . 1 , ^1 = ^2 = 70. lu Flgurc 2, we see that with small 
noise intensities the simulation of the two models are very similar; with higher small 
noise intensities, the simulations are very different. This is due to the fact that the 
behavior of the two diffusion processes near the boundary "6 = 0" are different: 
with the model (4) the washout regime is attainable which is not the case with 
the model (19). In Figure 3 we compare the evolution of the washout probability 
t — )■ F{Bt = 0) for both models, we clearly see that the model (19) does not give 
account for this probability. 

5.2 Simulation with the Haldane growth rate function 

In this test we use the Haldane growth rate function (3) and the parameters: k = 2, 
= 2.4 (mg/1), D = 0.1 (1/h), /i = 5 (1/h) , fc, = 10 (mg/1), a = 0.03: Ci = C2 = 
0.01. The initial law is (5o,5o) ~ A/'(1.5, 10"^) Ar(0.68, 10"^). The discretization 
parameters are s^^^ = 3, fe^ax = 2.5, 6 = 0.25, Ni = N2 = 300. 
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case l.a 



substrate 

case 2. a 



substrate 

case l.h 



substrate 

case 2.h 



Figure 2: In cases "1" the diffusion coefficients are cri(s) = ci a/s and (72(6) = C2 v^; 
in cases "2" the diffusion coefficients are cri(s) = ci s and ct2(6) = C2 &. In cases 
"a" Ci = C2 = 0.005; in cases "b" Ci = C2 = 0.02. For small noise intensities 
"1" and "2" behave rather similarly. For higher noise intensities 
(cases "b"), as the law -Kt of {St,Bt) is closer to the absorbing "washout" boundary 
{(s, b) G M^; = 0}, cases "1" and "2" behave rather similarly. See Figure 3 for the 
evaluation of the washout probability. 
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0.08 



■£ o.( 



Q- 0.( 



model 1 
model 2 




Figure 3: Washout probability 
for the case "1" (model 1: <Ji{s) 

(Ti(s) = CiS, (72(6) = C26). 



Following Figure 2: we compute t — )■ 
■ Ci ^/s, o"2(&) = C2 y/b) and for case "2" 



. = 0) 
model 2: 





1 1.5 

substrate 



Figure 4: Phase portraits for the system (1) for the Monod growth function (left) 
and the Haldane growth function. Left (Monod case): there are two equilibrium 
states: the washout equilibrium (red dot) is unattractive, the equilibrium point 
{s*,b*) with s* = fcs -D/(/i„,ax - D) (solution of /x(s) = D) and b* = (sj^ - s*)/k is 
attractive. We suppose that /i,„ax > D. The dashed line is 6 = (Si^ — s)/A;, in blue two 
trajectories (blue circles: initial positions). Right (Haldane case): the washout 
is still an equilibrium point but now it is attractive, there are two other equilibrium 
points given as solutions of fi{s) = D (we suppose that it admits two separate 
solutions), {s\,b\) is attractive (corresponding to the smallest value of s), (33,62) 
unattractive. The black dashed curve separates the two basins of attraction; in blue 
four trajectories (blue circles: initial positions). 
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Figure 5: Evolution of the distribution law of Xf. for each time t, the density pt{s, b) 
together with the washout density qt{s); the dashed curve separates the two basins 
of attraction. The mean of Xq is on this curve. See comments in the text. 
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In Figure 5 we plot the time evolution of the distribution law of Xf. for each time 
t, we represent (the approximation of) {pt{s, b); (s, b) G (0, s^^^) x (0, b^^^)) together 
with (the approximation of) {qt{s); s G (0, s,^^^)). In this test the mean of Xq is on 
this curve that separates the two basins of attraction (dashed white line): hence 
part of the mass will be attracted by (s*, 6^) and the other part will be attracted by 
the washout (sj^, 0) (see Figure 4). 

For t = we plot all the trajectory {x(t))teio;80] (white line). At the beginning 
the distribution law starts to "stretch" between the two attractors {t = 24). At 
t = 32, part of the mass is already on the point {sl,bl). Note that at this instant 
Pt{s,b) is bimodal and x(t) is a good approximation of K{Xt), but it is a poor 
statistics for Xt. At the final time t = 80, the deterministic trajectory x{t) reaches 
the equilibrium point {sl,bl) and 13% of the mass has been trapped by the washout 
absorbing boundary and some mass is still in the washout basin and will be trapped 
by the boundary "6 = 0" . 

A General finite difference scheme for n-dimensional 
diffusion processes 

Let Xt be the following diffusion process: 

dXt = b{Xt) dt + aiXt) dWt 

where Xt takes values in M", 6 : M" i-)- M", a : M" h-)- M"^™, and Wt is a standard 
Brownian motion with values in M™. Let a = a a* : ^ M"^". The coefficients 
are supposed to be locally Lipschitz and at most of linear growth. 

The probability density function p{t, x) of Xt is solution of the following Fokker- 
Planck equation: 

^^p{t,x)=C*p{t,x) (20) 
where C is the infinitesimal generator defined by: 

n 1 " 

i=l «J=1 

We consider finite difference schemes based on the following stencil (for the compo- 
nents (xj, Xj)): 
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-6 o o- 



-o- 



-o o- 



-o- 



We use the following up-wind scheme [8]: 



4>(x+hi ei)-4>{x) 
hi 

4i{x)-<t>(x-hi ej) 



if /.(a;) >0, 
if <0. 



{x+hi Ei+hj ej) — (j>(x+hi e^) (f){x+hj ej)—<j){x) 



{x) — (p{x—hj Bj) (^(x— hi ei) — (f)(x—hi Bi — hj ej) 



if aij{x) > , 



2 hi 



(x+hi ei)—4>{x+hi Si — hj ej) <p(x)—4>{x—hj Cj) 



hj hj 
[x+hj ej)—{j>{x) ^i+hj ej)—<j){x—hi e^) 



if aij{x) < - 



for i, j = 1, . . . ,n, i ^ j . The last non-diagonal second order schemes correspond to 
the following diagrams: 



O o- 

t— t 

-o • o- 

t— t' 

-o o 

a„ (,T) > 



o o 

t— t 

o • 

t— t 

— o o- 



<iij(i-) < 



With notation f~^{x) = max(/(a;),0) and / (x) = max(— /(x), 0), we get the 
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following approximation: 

i 

+ H ["^(^ + hi ^i) - 2 '^(^) + (t'ix - hi Si)] 

i 

+ ^ { i}-^^^ + hiei + hj Cj) - (p{x + hi e^] - + hj ej) - ^{x)] 

+ [(l>{x) - (l){x - hj ej)] - [(f){x - hi ei) - (f){x -hiCi- hj Cj)]^ 



(^[(pix + hi ei) - (j){x + hiei- hj ej)] - [(t){x) - (j){x - hj Cj)] 

+ [(f){x + hj ej) - (f){x)] - [(p{x - hiei + hj ej) - ^{x - hi Cj)] j | 

-9\x)<y ^. 2^.1 l^i^j-i^j 2hihj } 

+ Hx + e.) { ^ + - E,;,^. } + E, E.,^, Hx + hj ej) ^ 

+ ^{x - e.) { m + _ feif } + E, E.,^, 0(x - e,) 

+ Ei { i^i^ + hiei + hj ej) + (p{x - /ij - hj ej)] 

+ YH~T~ i^i^ + hiei- hj ej) + (j){x -hiei + hj ej)] | 
the symmetry Ojj = Qji leads to 

+ E. Hx + h,e.){il^ + ^-0^- E,,^. ^ } 

+ "Pl^; e^} I + l^rdi^i YKjii / 

+ Ejj;i^j { St^J [<^(^ + /ii ei + /ij Cj) + <p{x -hid- hj ej)] 

+ ['?^(^ + ei - hj ej) + ^{x -hiei + hj ej)] } 
We get the following approximation of the infinitesimal generator: 

C(p{x) ~ Ch(p{x) = Ch{x,y) (p{y) 



yeGh 
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for X & Gh where Gh — {x = (fci hi, . . . ,kn hn) ; ki — 0, . . . , Ni, i — 1, . . . ,n} and 

_sr^n \fi{x)\ v^n f au(x) \ 
l^i=l hi Z^i=l\ hf '^j¥=i 2hihj / ' 

Ch{x, x-hiCi- hj Cj) = for j , 

Ch{x, x-hiCi^ hj Cj) = 1^ for i^j, 
otherwise. 



Ch{x,x) 
jCh{x,x± hi ei) 
Ch{x, x + hiei + hj Cj) 
Ch{x,x + hi Ci - hjCj) 

^h{x,y) 



B Boundary conditions for the finite difference 
approximation 

For the boundary points Gh \ Gh of the grid, we use the foUowing schemes: 
• For X e {{s, b) eGh] s ^0,be {0, b^,^)} 

Ch{x,x) -- 



hi /i2 



hi 



hi ' 



Ch{x,x - hi ei) 
Ch{x,x±h2 62) 



hi 



2 hi 



because /i(0, 6) = D S;,, , (71(0, 6) = , 



fl(x) ^ al(x) 
h2 



+ ^ , note that /2(0, 6) = -£> 6 < , 



Ch{x,y) =0 otherwise. 

For a; e {(s, 6) e G/j ; s = s^^, 6 e (0, 6^ax)} 



^) ^ _LM£)1 _ \Mx)\ _ ^ _ ^ 



hi /i2 /jf 

£/i(a:, x + /ii ei) =0 (set artificially to 0) , 

Ch{x,x-hiei) =^ + 4r' 
Ch{x,x±h,e,) + 

Ch{x,y) =0 otherwise. 

For X e {(s, 6) e ; s e (0, s^J, 6 = 0} 



£/i(a;,a; ± /ii ei) 



/f (^) , '^f(^) 



2 /if ' 



Ch{x, X ± /i2 62) = ^ + = because ^(s, 0) = (72(5, 0) = , 
jCh{x,y) — otherwise. 
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For X e {{s,b) e Gh] s e (0, s^,,), b = 6^,,} 



h2 



X + /i2 62) =0 (set artificially to 0) , 
Ch{x,y) =0 otherwise. 



For X = (0, 0) 



Ch{x,x + /ii ei) 
Ch{x,x - hi ei) 

X ± /l2 62) 



hi 



= because /i(0, 0) = D s,^ , cti(0, 0) = 
= because /2(0, 0) = (72(0, 0) = , 



= 



otherwise. 



For X = (s^ax, 0) 



Ch{x,x + hi ei) 
Ch{x,x - hi ei) 
£h{x,x± h2 62) 



l/iWI '^fW 



hi hf ' 

(set artificially to 0) , 

1/1 Wl , ^'ii^) 



= because f2{s^^ 
— otherwise. 



^0) 



CT2[S„ 



,0) =0 



For a; = (0, b^,^) 



jCh(x,x) 

Ch{x,x + hiei) 
jChix,x - hiei) 
Ch{x,x + /l2 62) 
jChix,x - /l2 62) 



hi 

hi 



ft2 



= because /i(0, b^J) = -D Si„ , (7i(0, fo^^x) 
= (set artificially to 0) , 



\Mx)\ _^ alix) 



h2 



hj ' 



jCh{x,y) — otherwise. 
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For X (^max? ^max) 



Ch{x,x + hi ei) 
Ch{x,x - hi ei) 
Ch{x,x + /iaea) 
Ch{x,x - /i2 62) 



l/iWI I/2WI -^fw 



hi /i2 

(set artificially to 0) , 

1/1 Wl , dM. 

hi ~^ h'f ^ 

(set artificially to 0) , 



/12 



otherwise. 
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